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Abstract. We give computational results to study the accuracy of several quasicontinuum meth- 
ods for two benchmark problems — the stability of a Lomer dislocation pair under shear and the 
stability of a lattice to plastic slip under tensile loading. We find that our theoretical analysis of the 
accuracy near instabilities for one-dimensional model problems can successfully explain most of the 
computational results for these multi-dimensional benchmark problems. However, we also observe 
some clear discrepancies, which suggest the need for additional theoretical analysis and benchmark 
problems to more thoroughly understand the accuracy of quasicontinuum methods. 



1. Introduction 

Multiphysics model coupling has captured the excitement of the engineering research community 
for its potential to make possible the numerical simulation of heretofore computationally inaccessible 
multiscale problems. The capability to assess the accuracy of these multiphysics methods is crucial 
to both the verification of existing methods and the development of improved methods. 

During the past several years, a theoretical basis has been developed for estimating the error 
of atomistic-to-continuum couphng methods [2l[3l[8Hia[I3[ia[2lH271[29l[Ml[SH36]. However, the 
accurate computation of lattice instabilities such as dislocation formation and movement, crack 
propagation, and plastic slip are primary goals of atomistic-to-continuum coupling methods; and 
a theoretical analysis of the accuracy of atomistic-to-coupling methods up to the onset of lattice 
instability of the atomistic energy has thus far only been achieved for one-dimensional model prob- 
lems [5l llH [T2 t ll4 y 34 y 35j. This theoretical analysis and corresponding numerical experiments have 
given a precise understanding of the varying accuracy of these methods for simple model problems, 
but it is not known to what degree these errors are significant for the multi-dimensional problems 
of scientific and technological interest. 

We have developed a benchmark test code to study the atomistic-to-continuum coupling error 
for a face-centered cubic (FCC) crystal. Following the benchmark investigation of Miller and 
Tadmor [31], the displacement U (xi, X2, X3) of the atoms from their reference lattice positions is 
constrained by the atomistic analogue of continuum "plane strain" symmetry: 

U = {Ui, U2, 0) and U (xi, X2, X3) = u (xi, X2) 

where the crystal coordinates are given in terms of the basis vectors defined using Miller indices [Ij 
by ei =[110], 62 =[001], and 63 =[110]. We note that the computation of the atomistic energy 
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and forces for an atomistic displacement with plane strain symmetry requires the summation of 
the interaction of each atom with neighboring atoms in three-dimensional space. 

We have recomputed the benchmark study of Miller and Tadmor [31j for the Lomer dislocation 
dipole to more clearly separate the errors due to continuum modeling error, atomistic-to-continuum 
coupling error, and solution error. We allow a more general atomistic domain that is fully sur- 
rounded by a continuum region. We have also investigated plastic slip for tensile loading. Further 
investigation is underway to study the coarsening error in the continuum region and the far-field 
boundary condition error. 

We have found that the patch-test consistent quasi-nonlocal (QNL) method [30] gave a more 
accurate critical shear strain than the popular "ghost force correction" quasicontinuum coupling 
method [301138] tested in [31] (see Tabled and Figure [Tj). This result confirms our earlier theoretical 
analysis of the dependence of the critical shear strain error on the accuracy of the atomistic-to- 
continuum coupling [14]. 

The motivation for the ghost force correction method (equation ()2.12p below) was to correct the 
large coupling errors of the original quasicontinuum method (QCE) [H] by applying a "dead load" 
correction |30 1 I31 1 [38]. We have recast the ghost force correction method in a numerical analysis 
setting as an iterative method to solve the equilibrium equations for the force-based quasicontinuum 
(QCF) method [319], and we have proven for a one-dimensional model problem that the ghost force 
correction (QCE-QCF) method is inaccurate near lattice instabilities because it uses the patch- 
test inconsistent QCE method as a preconditioner for the QCF equilibrium equations [15]. The 
benchmark tests we present here confirm that the QCE-QCF method is not as accurate as the QNL 
method near lattice instabilities. 

By contrast, we found that the error was smaller for the QNL method than for the QCE- 
QCF method far from lattice instabilities (see Figures [3l9]) , which is contrary to expectations based 
on our analyses of a ID model in |lllll3t [T5 | l34j. where we have shown that the QCF method (the 
limit of the QCE-QCF iteration) is a more accurate approximation than the QNL method. 

Another result of our numerical experiments that contradicts our ID analysis [15j is that near a 
slip instability the QCE-QCF method has comparable accuracy to the QNL-QCF method, which 
uses the QNL energy as a preconditioner. To explain this, we note that our ID analysis in [15j can 
be considered a good model for cleavage fracture, but not for the slip instabilities studied in the 
present paper. We are currently attempting to develop a 2D benchmark test for cleavage fracture 
to demonstrate that the QCE-QCF method may be less accurate than the QNL-QCF method near 
some types of lattice instabilities. 

2. The atomistic and quasicontinuum models 

2.1. A model for plane strain in the face centered cubic lattice. Let £ be a face centered 
cubic (fee) lattice [1] with cube side length a and nearest-neighbor distance a/\/2 . Our plane 
strain model is most easily derived by viewing the fee lattice C as being generated by the primitive 
lattice vectors 

ai = ^(0, 1, 1), 02 = ^(1, 0, 1), and as = ^(1, 1, 0). (2.1) 
The cubic supercells are then generated by the cubic lattice vectors 

Ai = a(l, 0, 0) = -ai + 02 + as, A2 = a(0, 1, 0) = ai - a2 + as, 
As = a(0, 0, 1) = ai + a2 - 03. 

We let P be the orthogonal projection of £ onto the plane with normal given by ai— a2 = f (— 1, 1, 0) 
(or the (110) plane using Miller indices since we have that A1—A2 = — 2(ai— 02)), and let Q = I—P 
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Figure 1. Illustration of the FCC lattice and the projection onto the (1 1 0) plane. 
The grey shaded region represents the (1 1 0) plane, the asterisks are the points in 
the triangular lattice 7", and the disks are the atoms in the FCC lattice £. The 
vertical lines emanating from the points in T represent the columns of atoms in the 
FCC lattice which project onto points of T- The reader should imagine that there 
is a vertical column of atoms emanating from each point of T. 

be the projection onto the line parallel to ai — a2 = f (— 1, 1, 0) (or the [1 1 0] direction using Miller 
indices pj|). We observe that the projection of C onto the plane normal to ai — 02 = f (— 1) 1; 0) 
(or the (110) plane using Miller indices) is a triangular lattice T := PC. Each point in the 
triangular lattice T := PC is thus the projection of a column of atoms with spacing parallel 
to ai — 02 = f (— 1, 1) 0) (or the [1 1 0] direction using Miller indices) as depicted in Figured) 

Now let a; be a finite subset of the triangular lattice T, and let be the right cylinder over uj in 
the fee lattice C. We will call a displacement f7 : — t- periodic if 

U{X) = u[x + n{ai - 02)^ for ah X € 17 and aU n G Z. 

We note that the scale of the periodicity is |ai — 02] = a/^/2. Any periodic displacement U reduces 
to a displacement ti : w — )• M'^ defined by 

u{x) := U{X) for any X G with PX = x. 

We let U denote the space of displacements of w. 

We will now derive an energy per period : U ^ M. and a corresponding force field on cj. Let (p 
be a pair potential. For example, we use the Morse potential in our numerical experiments below. 
We will adopt a fourth-nearest neighbor pair-interaction model. For this choice of the interaction 
model the quasi-nonlocal coupling method [40j can be applied. We note that, throughout, we 
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understand n-th nearest neighbours in terms of Euclidean distance in the reference configuration 

in M3 (i.e., in the 3D FCC lattice) 00]. 

Let A4 denote the C-interaction range, by which we mean the set of all vectors pointing from an 

atom in £ to one of its first, second, third, or fourth nearest-neighbors. Representative first, second, 

third, or fourth nearest-neighbor vectors are given in terms of the primitive lattice vectors (12. ip 

and cubic lattice vectors (j2.2p by 

/ a a\ . , / a a\ 

ai = (^0, -, -j , ^i = (a, 0, 0), a2 + a3 = [a, -, -j , 2ai = [0, a, a) , (2.3) 

with the full set of first, second, third, or fourth nearest-neighbor vectors given by symmetry. 
For X gQ, we let 

Mx — {BGM-.X + Ben}. 

We will often call directions B G Mx bonds. We note that Ai ^ Mx when X lies in a shallow 
layer along the boundary of O. We also note that we have Mxi = Mx2 if P^i = PX2 = x, hence 
we can unambiguously use the notation Mx '■= Mx where PX = x for x €z uj. 

For any B £ M and any displacement U{X), we define the difference operator by 

DBUiX) := U{X + B)- UiX), 

whenever X, X + B G Q. We define the set of projections b = PB of the bonds in M and Mx 
onto the plane (110) by 

1Z:=PM and lZx:=PMx- 
For periodic displacements U{X), we define 

Di,u{x) := u{x + 6) - u{x) = U{X + B) - U{X) = DbU{X) 
where u{x) = U{X) for x = PX and b = PB. It will also be convenient to use the notation 

Dbu{x) := Dbu{x) = DbU{X) for PB = b. 
The energy per period : U ^ M is given by 

^'^W^=IE E H\B + DBn{x)\) 

where we note that the difference in the deformed positions of atoms at X + B and X is B+Dbu{x). 
The corresponding force at x G is then given by 

F'^{u){x) := -Q^f{n){x) = E '^d^ + Dbu{x)\), (2.4) 

where -q^t^ above denotes the partial derivative with respect to the displacement u{x) of the atom 
with projected reference position x G uj. We define the first variation (or, gradient) of £"'{u) by 
6£°'{u), which is given for each x G a; by 

6£'^{u){x) := ^^e-{u){x). (2.5) 

We can then express ()2.4p in operator notation by 

T^{u) = -58'' {u). 

We will now formulate a boundary value problem for the energy We let F C w denote the 
part of the "boundary" where the displacements are constrained. In the numerical experiments 



A COMPUTATIONAL AND THEORETICAL INVESTIGATION OF QUASICONTINUUM METHODS 



5 



described below, T will be a shallow layer of atoms along all or some of the boundary of uj. Let the 
space of admissible displacements that are equal to uq £U on the boundary F be denoted by 

Adm(uo) := {u £ U : u{x) = Uq{x) for all x G F}. 

We will study the minimization problem (or more precisely, the local minimization problem): 

Find u G arg local miniS"(f). (2-6) 

DfEAdm(Mo) 

The Euler-Lagrange equation corresponding to problem (j2.6p is 

-5£''{u){x) = for all X G w \ F, 

u{x) = uq{x) for all x G F. 

In our numerical experiments, we consider only initial displacements uq satisfying uq : a; — ?• 
(1 1 0), where we recall that (1 1 0) is the plane with normal given by oi — 02 = f (— 1, 1, 0). We 
will call any displacement n : w — )• (110) a plane displacement, and we will let V C denote the 
space of plane displacements. For plane displacements, £°' can be interpreted as the energy of the 
two dimensional crystal lo computed using a pair potential which depends on the bond direction. 

The restriction of the energy per period £"'{u) : Z^/ — )• M to plane displacements V can then be 
given by 

where (ph : (0, 00) — )• M is given by 

Mr):= J2 4>{{r' + \QB\'f') 

{BeMa:.PB=b} 

where we recall that Q = I — P is the projection onto the line parallel to oi — 02 = f (— 1, 1, 0). 
Roughly speaking, <j)ij{\b+ Dhu{x)\) is the energy of the interaction of atom x with all of its neighbors 
in the column over x + b. For x G w, we will call the inner sum 

■■=h E M\b + DMx)\) 

beTlx 

the atomistic energy at x. 

It will be convenient to establish coordinates adapted to the triangular lattice T, which lies in the 
plane defined by the normal oi — 02 = f (— 1, 1, 0) (or the (1 10) plane). This plane is spanned by 
the orthogonal vectors of length in the [1 1 0] direction and of length a in the [0 1] direction 

given by 

yi: = i(yli+^2)=a3 = f(l, 1, 0), 
V2 ■■ = A3 = ai + a2 - = a(0, 0, 1). 

Throughout the remainder of this paper, we will denote the coordinates of (xi, X2) in the {Vi, V2 } 
basis by 

{1^1, U2) ■■= viVi + V2V2. 
The triangular lattice T is then generated by the basis vectors (1, 0) and (1/2, 1/2). 
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It can be checked that there are four distinct symmetry-related projected interaction potentials, 
(pbir), corresponding to interactions with bonds B such that 

= (1, 0) PB = bi and QB = n{ai - ag) for n = -1, 0, 1, 

62 = (2, 0) PB = 62 and QB = 0, 

64 = (1/2, 1/2) PB = 64 and QB = n(ai - ag) for n = -3/2, -1/2, 1/2, 3/2, 

65 = (3/2, 1/2) PB = 65 and QB = n{ai - 02) for n = -1/2, 1/2. 

where the planar bonds bi are displayed in Figure [2j The remaining interaction potentials, (/>b(r) 
for b E TZ, can be obtained by symmetry from (pbiir) = (pbaii") = (j)bj{r), (jyb^ir), (t^biif)-, and (^^^(r). 
We note that the lattice constant a is contained in the definition given in the previous paragraph 
for the coordinates {i'i,i'2)- 

Now observe that for u G V we have 8£"'{u) G V where the gradient 68'^ {u) is defined by (j2.5p . 
Therefore, if a gradient-based minimization algorithm (such as the steepest descent method or the 
nonlinear conjugate gradient algorithm discussed below) is started from a plane displacement, it 
will terminate at a plane displacement. Thus, since we will only consider initial configurations 
which are plane displacements, we will be able to use the restriction of to V H Adm(uo) in our 
numerical experiments. This is important since restricted to V H Adm(no) is the energy of a 
two dimensional crystal, and in that case the patch-test consistent atomistic to continuum coupling 
method developed in |37J can be used. 

2.2. Energy-based quasicontinuum approximations. We will construct a local approximation 
to the energy Our methods follow [40j. First, we define an extrapolated difference operator D. 
The difference operator D will approximate the vector pointing from an atom to one of its neighbors 
using only the vectors pointing from an atom to its nearest neighbors. The nearest neighbors are 
the atoms depicted in Figure O Recall that IZ := PM\ we thus have bi £ TZ for the vectors bi, 
z = 1, . . . , 7, depicted in Figure O and moreover, TZ is obtained by applying the lattice symmetries 
of T to 61, . . . , 67. For these vectors, we define the operators 

Dfe, fori G {1,3,4}, 

2Db„ 



Db, 
Db, 
Db, 
Dbe 
Db, 



Db,+Db„ 
Dbs + Dbr^ , 
2Db,, 



and extend the definition by symmetries of T so that Db is defined for all b £ TZ. 
Using the diff'erence operator D, we define the Cauchy-Born energy by 

S^'iu) :=1^ ^ 0,(|6 + ^,^(x)|). 



We will call the inner sum 



the continuum energy at x, so 



:=i M\b + Dbuix)\) 



bell:. 



XdU) 

The following remark gives a more detailed motivation for these definitions. 
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Figure 2. In the above display of the (110) plane, the origin is marked with a 
dot, the nearest neighbors of the origin are marked with asterisks, and the other 
neighbors are marked with circles. The bond vector bi is the vector pointing from 
the origin to the atom marked bi. 



Remark 2.1 (The Cauchy-Born Approximation). We call £'^^{u) the Cauchy-Born energy since, 
if the displacement u is interpreted as a piecewise linear spline with respect to the canonical 
triangulation of T, it can be rewritten as an integral over a stored energy density. 

To make this precise, we observe that the continuum energy at x G cj such that x + b ^ lo for all 
6 G 7^ is given for uniform displacements u^{x) := Fx — x where F £ M^^-^ by 

ben 

We thus define the Cauchy-Born stored energy density, Wcfe : K^""^ ^ 1^ U {+00} , by 

W,b{F) 



where v is the area associated with each atom (in the 2D lattice). Moreover, for u £U,let u denote 
the piecewise affine interpolant of u with respect to the canonical triangulation of T, let Vn denote 
the resulting displacement gradient, and let O denote the union of those elements. 

For periodic boundary conditions, or by modifying the functional £'^^{u) at the boundary P of 
the domain a;, one can use Shapeev's bond density lemma [37] to prove that 



We note that this modification of the functional £'^^{u) at the boundary does not affect the minimiza- 
tion of quasicontinuum energies whose set of admissible displacements, Adm(tio), are constrained 
in the entire boundary of oj, such as for the Lomer Dislocation Dipole problem that we study. More 
detailed discussions and analyses of the Cauchy-Born approximation can be found in [^[TTlflS]. 
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Now let a; := AU C be a partition of w into an atomistic region A and a continuum region C. 
We define the energy based quasicontinuum (QCE) energy by 

f''-(^):=5]f,"(u) + 5]ff(n). (2.8) 

xeA x&c 

Following |40j . we also define the quasi-nonlocal (QNL) energy by 

£'i^\u) ■■=\Y.Y. Xb{x)Mb + D,u{x)\) + (1 - Xh{x)) Mb + DMx)\), (2.9) 

where 

, . J 1 if X £ C and x + b £ C, 

Xb[x) ■— <^ Q otherwise. 

The reason for the introduction of the QNL method was that the QCE method did not pass 
the patch test [11^138]. that is, the coupling mechanism defined in (12. 8p results in non-zero forces 
(the "ghost forces") at the atomistic/continuum interface under homogeneous displacements (or, 
deformations). By contrast, it was shown in [40] that the QNL energy ()2.9p does pass the patch 
test |lllll4ll40j . as long as two atoms interact only when they share a common nearest neighbor. 
The authors of |40j consider only the case where the pair potential does not depend on the bond 
direction. Nonetheless, the argument which they give may be applied to show that the QNL 
energy (12. 9p passes the patch test also in the present case. 

The importance of passing the patch test lies in the fact that the "ghost forces" can be understood 
as a consistency error, which results in an 0(1) relative error in the displacement gradient [IT]. 
Moreover, it was shown in [14], that they can result in a 0(1) relative error in the prediction of 
critical loads at which lattice instabilities occur. 

To achieve an efficient quasicontinuum method, the positions of atoms in the continuum region 
are normally further constrained by piecewise linear interpolation |30j . The development of an 
implementable and patch-test consistent quasicontinuum energy that allows coarsening by piece- 
wise linear interpolation in the continuum region has been achieved so far only for two-dimensional 
problems with pair potential interactions [37j. The development of patch-test consistent quasicon- 
tinuum energies for many-body potentials such as the popular embedded atom method (EAM) or 
for three-dimensional problems is an open problem. 



2.3. The force based quasicontinuum approximation. The force based quasicontinuum (QCF) 
approximation gives a patch-test consistent approximation for atomistic-to-continuum coupling 
even with coarsening in the continuum region [6l[8] . Force-based multiphysics coupling methods are 
generally popular because of their algorithmic simplicity. However, force-based coupling methods 
are known to give non-conservative force fields [Bill], and our recent research described at the end 
of this section has discovered additional stability problems for QCF approximations and solution 
methods |15]. 

Instead of approximating the energy £^ by a local energy, we can approximate the forces J-'{u) 
directly. This leads to the force based quasicontinuum (QCF) approximation. As above, let 
uj := ^ U C be a partition of oj into an atomistic region A and a continuum region C. Define 
jrqcf . y ^ y by 

iii^: (2.10, 
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We consider the boundary value problem 

7"'="=^ (u)(x) =0 for all X G w \ T, 

u{x) = Uo{x) for all x G F. 

It was shown in |12[|131 [28]. for one-dimensional model problems, that the solution of problem (j2.1ip 
approximates the solution of problem (j2.6p . It is also true that the QCF method passes the patch 
test [8j. 

In our numerical experiments we will solve equation (|2.11|) using an iterative method defined 
in [HllH]. Let uq G V be given, and let S'^" be either f"'"' or £1"". Then we let Un+i G Adm('Uo) n V 
be the solution of the problem: 

Mn+i £ arg local min 

^gAdm{«o)nV [ ^ ^ J 

The method (j2.12p with E^'^ = 8'^'^^ is the popular ghost force correction method proposed in [30] . 
We have proposed the method = S'^^^ as a ghost force correction method with improved stability 
properties pj- For = 8'>'"', we cah the iteration ([21^ the QCF-QCE iteration; for = 
we call it the QCF-QNL iteration. 

The Euler-Lagrange equation for problem (12.121) is 

58'^^{un+i){x) = F'^^f {un){x) + 58'^\un){x) for all X G w \ F, 
Un+i{x) = Uq{x) for all x G F. 

It is easy to see that, if the sequence of iterates converges, then its limit must be a solution 
of problem (j2.1ip . Moreover, under certain technical conditions related to the stability of the 
preconditioner 8'^'^, it can be shown that the sequence indeed converges |j8|[9lll5j. 

A QCF solution method that can theoretically give inaccurate results is the use of a nonlinear 
conjugate gradient solver for the force-based quasicontinuum equations [SIEIIEH]. We have proven 
that the linearization of the force-based equilibrium equations is not positive definite |13j . which 
implies that the conjugate gradient solution of this problem is unstable [15|. We have discovered 
from informal discussions with computational physicists that their conjugate gradient iterative 
solution of force-based multiphysics coupling methods sometimes oscillates rather than converges, 
a phenomenon partially explained by our theoretical analysis of this instability. We are developing 
benchmark tests to further study the reliability of the conjugate gradient solution of the force-based 
quasicontinuum equations. 

We have shown theoretically and computationally for one-dimensional model problems that the 
GMRES method is a reliable and efficient solver for the force-based quasicontinuum equations [15] . 
We are thus also preparing multi-dimensional benchmark tests to evaluate the reliability and effi- 
ciency of the GMRES solution of the force-based quasicontinuum equations. 

3. Numerical Experiments 

3.1. The Lomer Dislocation Dipole: Analytical Model. Following the numerical experiments 
presented in [31j . we consider a dipole of Lomer dislocations, as depicted in Figure [3l The Lomer 
dislocation has Burgers vector h = [110] in the (0 01) plane [20J. The dipole should be a stable 
equilibrium if the distance between the cores is large enough so that the Peach-Koehler elastic 
attraction of the dislocations is weaker than the critical force needed to overcome the Peierls energy 
barrier [31j . 
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7=applied shear 




[001] 



[110] 



L, 



Figure 3. Tlie atomistic and continuum (yellow) domains for the Lomer dislocation 
dipole. The size of the total computational domain is Li = 75xa/-v/2 and L2 = 60xa 
where a is the length of the sides of the cubic supercell. The atomistic domain is 
given by Di = k X a/\/2 and D2 = k x a for k = 3, 4, 5. The separation of the dipole 
used in the initial guess Ueias is w = 10 X a/V2. 



We will study the stability of the dipole when a shear is applied to the crystal. By a shear, we 
mean a homogeneous deformation of the lattice T that takes the form 

X ^ a{'~f)x where 17(7) := 

We call 7 the shear strain. We say that the shear is positive if 7 > and negative if 7 < 0. More 
generally, we will apply a shear to the deformation x + u{x) given by the displacement u{x) to 
obtain the sheared deformation 

x + u{x) I— (7(7) {x + u{x)) . 

A positive shear applied to the top and bottom boundaries of the crystal changes the energy land- 
scape to favor the movement of the dislocations apart. As 7 increases, the shear force overwhelms 
both the Peierls force and the Peach-Koehler elastic attraction. Thus, the original equilibrium 
configuration of the dipole becomes unstable, and the dislocations move away from each other. We 
call the value of 7 at which this instability occurs the critical strain. In our numerical experiments, 
we will simulate the process of slowly shearing the crystal until the critical strain is reached. 

To understand the movement of the Lomer dipole theoretically, we model the local minima of the 
atomistic energy for displacements constrained on the boundary T to be a{'j)x — x and constrained 
to be in the energy well of a Lomer dipole with separation w : 

_ inf S^iv) ^S^iw, 7) = Cisfit(li') + -^dipole attraction (W^) 

ii£Adm((T(7)z— w) 

~^ "-^applied shear T) ~^ "-^boundary effect 7) 

where Adm {a{^)x — x, w) is the set of admissible displacements roughly described above. Here 
the classical Peierls-Nabarro misfit energy [20j, the classical Peach-Koehler interaction energy |20] . 




A COMPUTATIONAL AND THEORETICAL INVESTIGATION OF QUASICONTINUUM METHODS 



11 



and our modeling of the effect of the appHed shear are given for w >0 and 7 > by 



fib 27TW 

gmisfitH = ^ .1 _ . cos-^exp 



ird 



{l-u)b 



(w)-i^ln- (3-1) 

'^dipole attraction V"' J — 27r L ' 

-^applied shear(^' T) = PlW + ^27 —, 

where is a characteristic shear modulus, is a characteristic Poisson's ratio, L is a characteristic 
length, b = a/\/2 is the magnitude of the Burger's vector, d = a/-v/2 is the interatomic spacing 
between the (110) planes, and /32 and (3^ are constants with /^s > 0. Although our model for the 
misfit energy '?misfit(^) ^^'^ dipole attraction <f^ipoie attraction ('"^) well-known approximations 
of the atomistic energy £°'{v) [20], our bilinear model for the shear energy '?3iear(''^' 7) derived 
from the atomistic energy and our justification is based only on the qualitative behavior it predicts 
for the stability of the Lomer dipole. 

After scaling the atomistic energy £°'{v), the dipole separation w, and the shear strain 7, and 
after neglecting the boundary effect "^boundary effect ~ obtain from (13. ip the model 

£"'{w,j)=cosw + (3ln'w — 'yw (3.2) 

for some /3 > 0. We chose 13 = 12 in our numerical experiments. We then find that the force 
conjugate to the dipole separation w is given by 

= sinzi; h7. (3.3) 

ow w 

We can see from the force field displayed in Figure H] that the dipole separation w becomes unstable 
at a critical shear strain 7 and will tend to infinity under gradient flow dynamics for the force 
field (|3.3p (for example, consider the stability of the equilibrium solution branch starting at w = 1.5 
as 7 is increased). The Lomer dislocations similarly separated to the boundary at a critical shear 
strain 7 in our numerical experiments (see Figure [5]). 

3.2. The Lomer Dislocation Dipole: Numerical Experiments. In out numerical experiments 
we use the Morse interatomic potential 

(p{r) := (1 — exp(a(r — 1)))^ for r > 0, 

with 

a := 4.4. 

We remark that we did not choose a = 4.4 in order to model some specific material. Rather, we 
chose this value of a as small as possible while retaining a stable dipole. For smaller values of a 
we were unable to compute a stable dipole, whereas, if a is chosen too large then the atomistic, 
Cauchy-Born, and quasicontinuum models are almost identical since the first neighbor interactions, 
which would become dominant, are treated identically in all models. 

We chose the lattice spacing a so that the FCC lattice £ is a ground state for the atomistic 
energy defined by the Morse potential and using a cut-off radius that includes first, second, third, 
and fourth-nearest neighbor interactions (determined from the reference positions). Specifically, we 
let a be the minimizer of 

V'(r) := 
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Force Field 




-5 1 1 ' ' 1 1 1 ' 1 

0.5 1 1.5 2 2.5 3 3.5 4 4.5 

w 



Figure 4. The force field ([33]) for the Lomer dipole model (132]) with /3 = 12. 

Separated Lomer Dipole Dislocations 




10 20 30 40 50 60 70 



Figure 5. Separation of the Lomer dislocations to the boundary after continuation 
past a critical strain to 7 = 0.070. Cartensian coordinates are used to label the axes. 
The top and bottom of the deformation of the reference domain uj are not displayed, 
but the entire width of the deformed reference domain is displayed. The triangular 
mesh was constructed from the reference lattice. 
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which is the energy per atom in the undeformed lattice with lattice spacing r since each atom in 
the £ lattice has 12 nearest neighbors, 6 next-nearest neighbors, 24 third-nearest neighbors, and 
12 fourth-nearest neighbors (see (|2.3p ). The minimum is attained when a = 1.3338. 

We let cj be the rectangular domain consisting of 150 columns of atoms parallel to V2 with each 
column containing 60 atoms defined in the (z^i, 1^2) coordinates by 



{{k, : < A: < 74, < / < 59} 

U {{k + 1/2, / + 1/2) : < < 74, < / < 59} . 

We take the boundary F of to be a layer of atoms four rows deep around the edges of cj, so that 
r = {{k, I) :k = 0,l, 73, 74, < / < 59} 

U {{k + 1/2, I + 1/2) : A; = 0, 1, 73, 74, < / < 59} 
U{(A: + l/2, ^ + 1/2) 
U{(A: + l/2, ^ + 1/2) 



< A; < 74, / 
< A; < 74, / 



0, 1, 58, 59} 
0, 1, 58, 59} . 



With this choice of F, for any atom x G w \ F, all fourth-neighbors of x belong to oj. Thus, there 
are no boundary effects for the reference lattice when a = 1.3338. This choice of F corresponds to 
taking Dirichlet boundary conditions on the entire boundary. Throughout the rest of this section, 
we will consider only displacements u so that for some 7 > 0, we have u{x) = a{'y)x — x for all 
X £ T. For such a displacement, we will call the shear strain 7 the shear on the boundary. 

We will now discuss how the stable Lomer dipoles displayed in Figure [6] are computed. We first 
construct an initial guess Ueias for the displacement field corresponding to a Lomer dipole by using 
the displacement fields of the isotropic, linear elastic solution for the edge dislocations [20] at the 
left and right pole of the dipole, respectively. More precisely, the formula for these displacement 
fields is 



'2^ 



tan 



_l X2 - X2 



(Xi - XI){X2 - X2*) 



Xl 



_^ _^ ^ ^ 

2{l-u){{xi-x\f + {x2-xlf^ 



U2 



2^ 



1 - 2v 
4(1 -z^) 



+ 



In ( Xl 



+ {X2 



*\2 



*\2 



(Xl - X^;)^ - (X2 - X2) 



4(l-z.)((xi 



+ {X2 



where * = L,R refers to the left or the right edge dislocation, respectively, b* is the corresponding 
Burgers vector, and x*, i = 1,2, are the components of the position vectors |20) . In our numerical 
experiments, we set the initial positions of the Lomer dislocations in the reference configuration to 
be (in the (1^1, 1^2) coordinates) 



{^1,^2) = ( 32,30 + 



and 



(-lV2^> 



42,30 



1 



We note that initial dislocation positions are placed between the center row of atoms at z/2 = 30 
and the first row above the center 1^2 = 30 + 1/2, with i/^" = 30 + 1/6 and z^l^ = 30 + 1/3 placed 
symmetrically about z^2 = 30 + 1/4 (see Figure [6|) . 



The displacements u ,u , are estimates of the displacement fields for isolated edge dislocations. 
By superimposing them, we obtain an estimate of the elastic displacement field for the dipole 
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(without shear): 

Uelas =U^ + U^. 

Finally, we apply a shear deformation to the deformation field x + Ueias , 

X + uliasix) = cr(7) (2; + Uelas{x)) , 

which yields an initial guess for the displacement field for the Lomer dipole under applied shear: 

Uelasix) = 0"(7) {X + Uelas{x)) - X. 

We use the displacement u^^^^ = cr(7o)^ie/as(2;) — 2; as a starting guess for a preconditioned nonlin- 
ear conjugate gradient method (P-nCG), which is described in Section [5l to solve the minimization 
problem (|2.6p . 

It is challenging to find a stable dipole. In fact, we were unable to find a stable dipole without 
applying a positive shear to the boundary. That is, we did not find a Lomer dipole that was a local 
minimum of the atomistic energy (j2.7p subject to the boundary condition u{x) = for all x G F. 
We were also unable to find a stable dipole when the parameter a in the Morse potential was too 
small. Essentially, we found that as a decreased, the interval of shear strains for which a given 
dipole was stable became smaller. After some experimentation, we did find a stable Lomer dipole 
with Q = 4.4 and with shear on the boundary 70 = 0.0375. 

Our atomistic and quasicontinuum models utilize a fourth-nearest neighbor cutoff calculated 
from the positions in the reference lattice, which is in fact the largest possible neighborhood for the 
QNL method applied to an FCC lattice [IQ]. Such a cutoff is acceptable for atomistic configurations 
that are "close" to the reference lattice, however, for the highly deformed positions near the dipole 
a cutoff in deformed coordinates ought to be used. Even though the dominant nearest-neighbour 
interactions are always included in our calculation, a more precise understanding of the error 
committed for second- to fourth-nearest neighbours is required. 

Next we describe how we simulate the shearing of the dipole. For simplicity, we will explain how 
this is done for energy-based methods before we discuss the force-based method. Let uq be the 
displacement field of the stable Lomer dipole discussed above. Suppose that 70 is the shear on the 
boundary for uq. Let 7^ be an increasing sequence of shear strains, and let d'jk+i ■= Ik+i — Ik- S, 
can be either the atomistic energy (12. 7p . the QCE energy (12. 8|) . or the QNL energy (12. 9|) . 

We perform the following iteration starting with the displacement uq. Suppose we have a dipole 
with displacement field which solves the boundary value problem 

— (5iS(nfc)(x) =0 for all X G a; \ F, 

"u-kix) = (T(7fc)x — X for all 2; G F. 

We let 

ui+i{x) = cr{hk) (x + Ukix)) - x. 
Observe that u^+i satisfies the boundary condition 

Uk+i{x) = a{jk+i)x — X for all x G F. 

We use the P-nCG algorithm to compute a new local minimizer Uk+i "near" u^j_i, which satisfies 
the new boundary condition Uk+i{x) = a{'yk+i)x — x for all x G F. We call this process the shear 
loading continuation. 

We terminate the shear loading continuation when we can no longer find an equilibrium Uk+i{x) 
close to the initial guess a{5jk+i) {x + Uk{x)) — x. In practice, this means that we stop when the 
nCG algorithm returns a deformation in which the cores of the dislocations have moved apart. If 
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the dislocations move apart at the {k + 1)*'' step of the iteration, then we assume that the critical 
strain must be between 7^ and 7^+1- 

For the force based approximation (|2.10[) . we use a similar iteration to find the critical strain. 
At each step in the iteration, we solve the boundary value problem 

= for all X e \ T, 

Uk+i{x) = a{'yk+i)x — X for all x G F. 

We solve equation ()3.4p using the ghost force correction iteration defined in equation ()2.12|) . We 
use the same initial guess 

ul+i{^) '■= cr((57fc+i) {x + Uk{x)) - X 

as in the energy-based case, and we solve each iteration of ghost force correction using the P-nCG 
algorithm. 

We presented computational results for the deformation of a Lennard-Jones chain under tension 
in [9] that demonstrate the necessity of using a sufficiently small parameter step size to ensure that 
the computed solution remains in the domain of convergence of the ghost force correction iteration 
(defined below as QCE-QCF) method. These results exhibit fracture before the actual load limit 
if the parameter step size is too large. We thus conclude that the shear strain must be 

sufficiently small to ensure that our computed solution remains in the domain of convergence of 
the QCE-QCF method since it would otherwise predict instability for the Lomer dipole before that 
predicted by the QCF method. 

We perform the shear loading iteration for each of the quasicontinuum methods using three 
different atomistic regions, as depicted in Figure [6l The atomistic region (3) (resp., (4), (5)) is a 
box containing the dipole such that there are three (resp., four, five) rows of atoms between the 
continuum region and each of the atoms in the pentagons surrounding the cores of the dislocations, 
or more precisely, the atomistic region [k) is the box given in {vi, V2) coordinates by [32 — /c, 43 + 
k] X [30 — A;, 30 + k] for /c = 3, 4, 5. Throughout the remainder of this section, QNL-QCF (resp. 
QCE-QCF) will refer to the QCF method implemented using the ghost force correction iteration 
preconditioned by the QNL (resp. QCE) energy [15]. We wih let QNL-QCF(k) (resp. QCE- 
QCF(k), QNL(k), QCE(k)) refer to the QNL-QCF (resp. QCE-QCF, QNL, QCE) method with 
atomistic region (k) for k one of 3, 4, or 5. 

We note that the atomistic regions studied in the benchmark tests given in [31| extended to the 
lateral boundaries of the computational domain w, or, following the notation used in Figure [3] these 
tests study the case 2Di + w = Li. Thus, the benchmark tests |31J do not present the capture of 
a Lomer dipole in an atomistic region fully surrounded by a continuuum region and do not study 
the effect of coupling error in the dipole plane (110). 

In Table [H we give the critical strains for each method in decreasing order. The columns 7" and 
7^ are the shear strains at the last step before the dislocations moved apart, and at the first step 
for which the dislocations moved, respectively. The column labeled "Error" displays the relative 
error from the critical strain of the atomistic model, calculated by the formula 

Error = I^^L- — Jsi^ (3.5) 

where ^at is the mean of 7" and 7+ for the atomistic method, and 7qc is the mean of 7" and 
7^^ for the quasicontinuum method. We wish to stress that if 7"*" for some method is equal to 7" 
for another method, then no strong statement can be made comparing critical strains of the two 
methods. For example, our data do not show that QNL(5) has a significantly higher critical strain 
than QCE-QCF(5). Rather, our data suggest only that the critical strain of the QCE-QCF(5) 
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(a) Atomistic region (3) (b) Atomistic region (4) 
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(c) Atomistic region (5) 

Figure 6. The atomistic regions for the Lomer dislocation dipole. Not all of the 
continuum region is displayed. The mesh is the Delaunay triangulation [7J of the 
reference lattice (i.e., with the atoms as the nodes). The atomistic region is shaded 
grey. 



method is in the interval [0.038172,0.038177] and that the critical strain of the QNL(5) method is 
in [0.038177,0.038181]. Thus, the critical strains of the two methods are very close. 

The reader will observe that the critical strain of the fully atomistic model is the highest. The 
critical strain predicted by the QCE method is the lowest, and is the farthest from the critical 
strain of the atomistic model. This is not surprising in light of the results on the stability of one 
dimensional chains obtained in |14j . It is surprising, however, that the QNL method predicts the 
critical strain of the atomistic model more accurately than the QNL-QCF force-based method, 
since we expected that the stability of the QNL-QCF method would be determined by the QNL 
preconditioner [15]. It also somewhat surprising that the QCE-QCF iteration predicted the critical 
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Table 1 . The critical strain for the Lomer dislocation dipole under shear. The Error 
is given by {'^at ~ Iqc) llat in percent which is defined precisely in the paragraph 
surrounding p.Sh . 



Method 


Critical Strain 




Error 


Atomistic 


7 

0.038190 


7+ 
0.038194 


0, 


,000 


% 


QNL(5) 


0.038177 


0.038181 


0, 


,034 


% 


QCE-QCF(5) 


0.038172 


0.038177 


0, 


,046 


% 


QNL-QCF(5) 


0.038172 


0.038177 


0, 


,046 


% 


QNL(4) 


0.038164 


0.038168 


0, 


,069 


% 


QCE-QCF(4) 


0.038155 


0.038159 


0, 


,092 


% 


QNL-QCF(4) 


0.038155 


0.038159 


0, 


,092 


% 


QNL(3) 


0.038111 


0.038116 


0, 


,206 


% 


QCE-QCF(3) 


0.038081 


0.038085 


0, 


,286 


% 


QNL-QCF(3) 


0.038081 


0.038085 


0, 


,286 


% 


QCE(3) 


0.037750 


0.037775 


1, 


,125 


% 



strain as accurately as the QNL-QCF iteration. The analysis given in [15j suggests that the QCE- 
QCF iteration should lose stability at a lower strain than the QNL-QCF iteration. 

Figure [7] consists of graphs showing the relative w^'°^ error versus the shear strain for various 
quasicontinuum methods. The w^'^ norm is defined by 

II II |-Dfen(x)| 
I 1^^,1,00 := max max , 

and we define the relative error in w^'°^ by 

/ \ I |Uqc ~ Ua| |wl'°° (oa\ 
errrel(Uqc,Ua) := — \ — . (3.6) 

Here Uqc denotes a solution of one of the quasicontinuum models, and Ua denotes a solution of the 
atomistic model with the same boundary conditions. 

The reader will observe that the error of the QCE method is the greatest. This is largely due 
to the presence of ghost forces [TU]. The error of the QNL method is the least. Again, this is 
surprising since various analyses [TT1[T31[TS1[251[31| suggest that the force based method should be 
more accurate than the QNL method. However, we note that a similar effect was observed in one- 
dimensional numerical experiments in [28]: while for large atomistic regions the QCF method was 
considerably more accurate than the QNL method, for smaller atomistic regions the QNL method 
was clearly more accurate. The atomistic regions studied here in our numerical experiments for the 
Lomer dislocation dipole are very small. 

3.3. Tensile Loading. In our next experiment, we examine the accuracy of various quasicontin- 
uum methods for the computation of the resistance to slip under an applied tension. By a tension, 
we mean a homogeneous deformation of the lattice T which takes the form 



X I— )• r(7)x where t(7) := 



1+7 
1 



when expressed in the basis of the coordinate vectors Y\ and V2. We will let t(7) denote such a 
tension, and we will call 7 the tensile strain. We call a tension positive if 7 > and negative if 
7 < 0. 
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Figure 7. The relative error errrei(uqc, Ua) 
the Lomer dislocation dipole. 
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In our numerical experiments, we use the same pair potential and lattice spacing as for the Lomer 
dipole problem. We let uj be the rectangular domain consisting of 120 columns of atoms parallel to 
V2 with each column containing 15 atoms defined in the (z/i, z/2 ) coordinates by 

uj = {(/c, /) : < A: < 59, < / < 14} 

U {{k + 1/2, / + 1/2) : < /c < 59, < / < 14} . 

We will now let V consist of four columns of atoms parallel to V2 at the left edge of cj, and four 
columns at the right edge of cj, that is, 

r = {(/c, /) : A; = 0, 1, 58, 59, < Z < 14} 

U {{k + 1/2, I + 1/2) : A; = 0, 1, 58, 59, < / < 14} . 

This corresponds to choosing Dirichlet boundary conditions on the sides of uj parallel to 1^, and 
free boundary conditions on the sides of oj parallel to Vi. 

To find the critical strain for the crystal under tension, we perform an iteration similar to the 
shear loading iteration. Let 7^ be an increasing sequence of tensile strains with 70 = 0. We let uq 
be the displacement field of the undeformed reference lattice oj described above. Then we let n^+i 
solve the problem 

F{uk+i){x) =0 for all X G w \ F, 

Uk+i{x) = r(7fc+i)x — X for all x G F. 

Here F is the gradient of one of the energies (|2.7p , (j2.8p , (|2.9p , or else the force-based method (|2.10p . 
We solve equation (j3.7p using the same methods discussed above for the shear loading iteration, 
except that in this case we start each step with the initial guess 



n^i(x) :=r(i±gf) {x + u,{x))-x. 



We call this iteration the tensile loading continuation. We stop the continuation when the P-nCG 
algorithm (or a ghost force correction iteration) returns a deformation in which a slip has occurred. 

Our numerical experiments are designed to test how well the boundary between the atomistic and 
continuum regions resists slip under tensile loading. When we perform the tensile loading iteration 
for the fully atomistic model, we found that a slip tends to occur along the slip plane indicated 



in Figure 8(a) when a critical tensile load is reached. Figure 8(a) shows the configuration of the 
crystal immediately after a typical slip has occurred. We note that the slip allows the crystal to 
accommodate an increased tensile strain with a decrease in the energy. We then chose an atomistic 
region A whose boundary is along that line, as depicted in Figure |8(b)[ 

We observe that the critical strains of the quasicontinuum models were lower than the critical 
strain of the atomistic model. Our results are summarized in Table [2j In Table [21 the columns 7" 
and 7"'' are the tensile strains at the last step before slip is observed and at the first step for which 
slip is observed, respectively. The column labeled "Error" is the percent error from the critical 
strain of the atomistic model. 

We observe that the critical strain of the QCE model is lower than the critical strain of any 
of the other models. This is in agreement with the one dimensional stability results established 
in |14] . However, we were surprised to find that the critical strain predicted by the QNL method is 
higher than the critical strain predicted by the QNL-QCF force based method, since we expected 
that the stability of the QNL-QCF method would be determined by the QNL preconditioner [15] . 
We were also surprised to observe that the critical strain of the QNL-QCF method is less than the 
critical strain of QCE-QCF. The one-dimensional analyses in [12[I15] suggest that the force based 
methods should have a comparable critical strain as QNL, and that the QCE-QCF iteration should 
lose stability at a lower strain than the QNL-QCF iteration. 
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(a) The natural slip plane. We choose an atomistic region whose boundary is 
the plane highlighted in grey. The mesh is the Delaunay triangulation ,7 of the 
reference lattice. 




(b) The atomistic region A is shaded grey. All atoms in A are marked with a blue 
dot. 



Figure 8. The atomistic {A) and continuum (C) regions for the tensile loading experiment. 



Table 2. The critical strain under tension. The Error is given by {'fat — Iqc) Mat 
in percent, which is defined in (j3.5p . 



Method 


Critical Strain 


Error 




7 


7+ 




Atomistic 


0.081635 


0.081649 


0.000 % 


QNL 


0.081593 


0.081607 


0.052 % 


QCE-QCF 


0.081242 


0.081270 


0.473 % 


QNL-QCF 


0.081101 


0.081130 


0.645 % 


QCE 


0.063900 


0.064200 


21.548 % 



Figure [9] consists of graphs showing the relative w^'°° error versus the tensile strain for various 
quasicontinuum methods. We observe that the QCE model had the greatest error. This is primarily 
the result of the ghost forces which arise in the QCE model. The QNL model had the lowest error. 
Again, this is somewhat surprising in light of the one dimensional theory |li y i3 1IT5ll281 l34j. which 
predicts that the force based methods should have lower error than the QNL method. 
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Figure 9. The relative error erri.ei(uqc, Ua) = ||uqc — Ua||wi>°°/||ua||wi.°° w^''^ for 
tensile loading. 
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4. Conclusion 

Materials scientists and engineers typically attempt to verify their multiphysics methods by 
benchmark tests. However, definitive conclusions from these benchmark tests are often not possible 
since they combine many sources of error (modeling error, coupling error, solution error, boundary 
condition error, etc.). Further supporting theory is needed to select a set of test problems that 
thoroughly samples the solution space. 

We are developing an improved theoretical basis for benchmarking atomistic-to-continuum cou- 
pling methods based on the multiscale numerical analysis developed by us and others. Our goal is 
to be able to reliably predict the accuracy of atomistic-to-continuum coupling methods for general 
deformations and loads from numerical experiments for a small set of mechanics problems. This 
set of mechanics problems should sample the fundamental modes of material instability such as 
dislocation formation, slip, and fracture. 

Our discussion of the benchmark tests presented in this paper give many examples of the pre- 
dictive success of the theoretical analysis we have developed during the past few years, but we 
also describe several cases where our theoretical analysis seems to predict a different outcome than 
our computational experiments. This discrepancy between theory and computational experiment 
occurs when our theoretical analysis does not adequately model the computational problem and is 
motivation to develop more general theoretical analysis. 

5. Appendix: A Preconditioned Conjugate Gradient Algorithm 

In this appendix, we describe the preconditioned nonlinear conjugate gradient optimization al- 
gorithm (P-nCG) and the corresponding linesearch method, which we used in the numerical exper- 
iments in this paper. 

Let £ : — ?• R U {+00} be continuously differentiable in {u : £{u) < 00}. If £ is an energy of 
the type discussed above, then the standard nonlinear conjugate gradient method [331 Sec. 5.2] is 
convergent, but very inefficient, due to the poor conditioning of the Hessian matrix at local minima. 

Let P G W^^^ be a symmetric positive definite matrix, e.g., a discrete Laplacian or a modification 
thereof. A simple but considerably more efficient method is obtained if all inner products in the 
conjugate gradient algorithm are replaced by a P-inner product {u,v)p = u^Pv, and all gradients 
6£{u) by P-gradients 6p£{u) = P~^5£{u). The P-gradient 6p£{u) "represents" the gradient 5£{u) 
in the P-inner product {v, w)p := v^Vw, since 

{5p£{u),w)p = {F~^5£{u))^Vw = 5£{ufw = {5£{u),w). 

In practice, we allow a new preconditioner to be computed at each step. A basic preconditioned 
conjugate gradient algorithm of Polak-Ribiere type can be described as follows: 

(0) Input: uo G M^; 

(1) Evaluate Pq; Qq = Vq'^5£{uq)] sq = 0; 

(2) Forn = 1,2,... do: 

(3) Evaluate P„; 

(4) On = Pn^5£{Un); 

(5) I3n = max{0, - c/n-i)p„/(9n-i,5'n-i)p„_i}; 

(6) Sn = -gn + l^nSn-i; 

(7) an ^ LINESEARCH; 

(8) Un = Un-1 + OnSn, 

In the following we specify further crucial details of our implementation: 
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(1) Preconditioner: The choice of preconditioner has the biggest influence on the efficiency of 
the optimization. We project all atoms onto a single plane and triangulate the resulting 
set of nodes. On this triangulation, we assemble the standard Pl-finite element stiffness 
matrix, K„, discretizing the negative Laplace operator. The preconditioner P„ is obtained 
by imposing homogeneous Dirichlet boundary conditions on the clamped nodes. 

(2) LINESEARCH: Our linesearch is implemented following Algorithms 3.5 and 3.6 in [33] 
closely. We use cubic interpolation from function and gradient values in the "Interpolate" 
step of Algorithm 3.6. We guarantee the strong Wolfe conditions, [33, Eq. (3.7)] with 
constants ci = 10~^, C2 = 1/2. 

(3) Initial guess for an: If the initial guess ai^"* for the new steplength an, which is passed to 
the LINESEARCH routine, is chosen well, then actual linesearch can be mostly avoided, 
which can result in considerable performance gains. Following [33^ Eq. (3.60)], and extensive 

experimentation with alternative options, we choose an^ = 2{£{un~i)—£{un~2))/{9n, Sn)p„- 

(4) Termination Criteria: We terminate the iteration successfully if the following condition is 
satisfied: 

(^||u„ - n„_i||oo < T0L2° or ||n„ - u„_i||p„ < TOL^^ 
and (||5£:(7x„)||oo <TOL^ or ||5„||p„ < TOLf) 

and {£{un-i) - £{un) < TOL^^ , 

where || • ||oo denotes the ^°°-norm, and || • ||p denotes norm associated with the P-inner 
product. The tolerance parameters are adjusted for each problem. Typical choices are 
TOL~ = TOL^ = 10-^ TOL~ = TOL^ = 10-^ TOL^ = 10"^. 

We terminate the iteration unsuccessfully if a maximum number of iterations is reached, 
or if the LINESEARCH routine is unable to make any progress. Note also that, since 
P-nCG is a descent method, we have £{un~i) — £{un) = \£{un-i) — £{un)\- 

(5) Robustness Checks: In addition, our algorithm uses various minor modifications to increase 
its robustness. Since these do not significantly affect its performance, we have decided to 
not give any details. 

The algorithm described above is both efficient and robust for most of the problems we consider. 
However, we wish to stress certain difficulties that arise in the presence of "meta-stable" states 
and particularly shallow local minimizers, which are difficult to distinguish numerically. In our 
experience, dislocations fall precisely into this category. 

On some occasions, our algorithm would fail when a particularly low tolerance setting was used. 
The reason for the failure is usually that the directional derivative along the search direction is 
non-negative (up to numerical precision) and hence the linesearch fails or stagnates. Replacing 
the conjugate direction by a steepest descent direction (one of our robustness checks) resolves this 
problem only partially. We were able to overcome these difficulties mostly by tweaking the various 
optimization parameters. 
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